Critical properties of the reaction - diffusion model 2A — > 3A, 2A 
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The steady-state phase diagram of the one-dimensional reaction-diffusion model 2A — > 3A, 
2A — > is studied through the non-hermitian density matrix renormalization group. In the absence 
of single-particle diffusion the model reduces to the pair-contact process, which has a phase transition 
in the universality class of Directed Percolation (DP) and an infinite number of absorbing steady 
states. When single-particle diffusion is added, the number of absorbing steady states is reduced to 
two and the model does not show DP critical behaviour anymore. The exponents 8 — v^/v± and 
/3/u± are calculated numerically. The value of P/v± is close to the value of the Parity Conserving 
universality class, in spite of the absence of local conservation laws. 
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I. INTRODUCTION 



Reaction-diffusion systems have attracted considerable 
interest in the past few years . While at sufficiently 
high spatial dimensions their critical behaviour is cor- 
rectly described by mean field rate equations, at dimen- 
sions below the upper critical dimension, where the effect 
of fluctuations becomes important this approach is not 
valid anymore. In this case one has to resort to other 
methods as, for instance, field theory Qj, exact calcu- 
lations via Bethe ansatz ||], Monte Carlo and cellular 
automaton simulations (see |2||| and references therein) 
or exact diagonalization techniques . Recently other 
types of approaches have been proposed as the density 
matrix renormalization group [9 and the standard real- 
space renormalization group |l0( . 

In this paper we study the critical properties of a one- 
dimensional reaction-diffusion model where the local dy- 
namics is given by the following rules. Consider a single 
species of particles (A) on an one-dimensional lattice. 
Each lattice site can be either occupied by a single parti- 
cle or be empty. A single particle may hop to an empty 
neighbouring site (diffusion). Two particles on neigh- 
bouring sites can either annihilate (2 A — > 0) or give birth 
to a new particle (2A — > 3A) in the case that one of the 
sites next to the couple is empty. The reaction rates for 
these processes are defined in equations (0J|^) below. 

Recently, Howard and Tauber |nj discussed a general- 
ized version of the above model where n particles may be 
created through the reaction 2A — ► (n + 2) A, and with 
an arbitrary number of particles per site. They employed 
a field-theoretical approach and found that the theory is 
not renormalizable. They argued, however, that for all 
n the model should not be in the directed percolation 
(DP) universality class. Their conjecture is based on the 
analysis of the associated master equation and on the 



massless nature of the field theory which is at odds with 
a phase transition of DP type. However, the critical be- 
haviour of the system, i.e. its universality class, remained 
unknown. In absence of diffusion, with at most one par- 
ticle per site, the n = I case is the Pair Contact Process 
(PCP), where the steady-state phase transition belongs 
to the DP universality class [|l2|,|lJ 



Here we present a non-perturbative study of the model 
corresponding to the n = 1 case of |rj| : the PCP model 
with additional single-particle diffusion, as obtained from 
density matrix renormalization group (dmrg) calcula- 
tions. The DMRG was introduced by White Q in 1992 to 
investigate numerically ground state properties of quan- 
tum spin chains. Because of its great accuracy, reliability 
and the possibility of treating large systems with a lim- 
ited computational effort, the dmrg has since been ap- 
plied to an ever increasing set of problems, as reviewed in 
p"5| . In particular, some recent studies were devoted to 
the application of the dmrg to non-hermitian problems 
1^,^-20 1 which appear frequently in many domains of 
physics, for example in low-temperature thermodynam- 
ics of spin chains and ladders, for models of the non- 
integer quantum Hall effect and in one-dimensional non- 
equilibrium systems. Reaction-diffusion systems belong 
to the latter class of models. New insight should be ex- 
pected from the application of the DMRG to them || . 

The paper is organized as follows: in Section |n| we 
define the model and recall some facts about reaction- 
diffusion systems. In Section III wc introduce a shift 



strategy to project out a trivial ground state from the 
quantum Hamiltonian. This improves the convergence 
of the DMRG. The method works if an exact expression 
for the eigenstate is available. In Section |^ we discuss 
the calculation of the critical points and exponents for 
the model and show that the model does not belong to 
the DP universality class. Section [V| concludes the paper. 
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II. THE MODEL 

We consider a one-dimensional lattice of length L with 
open boundary conditions, as usual in DMRG calculations 
Ea . The reactions occur with the following rates 
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and are parametrized by the diffusion constant d and the 
pair annihilation rate p. This defines the Pair contact 
process with diffusion (PCPD) model. 

For a first qualitative overview, we consider the mean 
field kinetic equations. To discuss the effects of diffusion 
we need the kinetic equations in the pair approximation, 
see e.g. [||. If n = n(t) is the spatially averaged single 
particle density and c = c(t) the spatially averaged pair 
density, we find 

n = -2(1 - d)pc+ (1 - d)(l -p) (n~c)c/n (4) 

2c + n (n-c)(c-n 2 ) 

c= -(1 - d)pc 2d 

n n(L — n) 

+(1 - d)(l - p) (n - c)(l - c)c/(n(l - n)) (5) 

In the d — ► I limit, non-trivial steady states only occur 
for c(t) — ► n(i) 2 and one recovers the single-site kinetic 
equation (rescaling time by a factor (I — d)) 



(I — p) n 2 (I — n) — 2p n 2 



(6) 



which expresses that particles may only be created on 
empty sites. The time-dependence of n(t) for t large is 



ioo + aexp(— t/r) 
l/(3p-I)-i- 1 



if P < Pc,mf(1) 
if p = p c , MF (l) 
if p > p c ,mf(1) 



(7) 



where ««, = (1 - 3p)/(l - p), r = (1 - p)/(l - 3p) 2 , 
Pc,mp(1) = 1/3 and a is a constant which depends on the 
initial conditions. At small p, where the creation process 
(|l|) dominates, the system is in the active phase with a 
non-vanishing particle density in the steady state. On 
the other hand, for p large, the pair annihilation process 
(^) dominates, the steady state particle density is zero 
and the system is in the inactive phase. In the entire 
inactive phase, and not only at the critical point, the ap- 
proach towards the steady state is algebraic, rather than 
exponential as found in most systems. 

The same kind of result also holds in the pair ap- 
proximation with d 7^ 1. In the steady state, we have 
c(oo) = n(oo)(l — 3p)/(l — p). The particle density van- 
ishes along the curve 
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< d < 1/7 
1/7 < d < 1 



(8) 



Close to the criticality, the particle density n(oo) ~ 
Pc,MF(d) — p, which is the same found from in the 
d — > 1 limit above. However, the pair density 



c(oo) 



Pc, MF (d) -p 
(PcM F (d) -p) 2 



d < 1/7 
d > 1/7 



(9) 



There should thus be different universality classes on 
both sides of the meeting ( "tricritical" ) point pt = 1/3, 
dt = 1/7, where n(oo) ~ (p t — p) 2 and c(oo) ~ (p t — p) 3 
imply modified critical exponents. 

The leading long-time behaviour along the critical line 
p = p c (d) is as follows 



i{t)~t- 1 / 2 , c{t)~t- 1 



if d > d t 



n(t) 
n(t) 



-l 
-l 



, c(i)~i~ 3 / 2 ;ifd = dt = l/7 (10) 
, c(t) ~ t^ 1 ; if d < d t 



while ~ t~ x and c(t) ~ t~ 2 in the entire inactive 
phase p > p c (d). In the active phase, the steady state is 
approached exponentially. 

Although generally believed to be qualitatively correct 
in sufficiently high dimensions kinetic equations 

cannot provide correct values of the exponents (and often 
not even the order of the transition) in low dimensions. 
In one dimension, the exact decay in the inactive phase 
is n(t) ~ i -1 / 2 (IT), where the exponent 1/2 is that of 
the process 2 A -> model ||, which is recovered tak- 
ing the limit p — ► 1 in the rates (0) and (||). However, 
the pair approximation suggests the presence of distinct 
universality classes for d large and d small. In addition, 
the exponent of the time-dependence of n(t) for d large 
equals 1/2 as found exactly for one-dimensional diffu- 
sion. This suggests that the upper critical dimension d* 
for that transition should be unity, see p^j , but the value 
of d* in the PCPD is not yet known. To what extent are 
the predictions of the pair approximation, in particular 
the existence of several distinct universality classes along 
the transition line, borne out ? 

An active state with a finite density of particles can be 
maintained only in the limit L — > oo. On a finite lattice 
any configuration of particles will decay towards an ab- 
sorbing configuration in a finite time. There are two pos- 
sible absorbing configurations which the system cannot 
leave: (i) the empty lattice and (ii) a lattice occupied by 
one single diffusing particle. If we take d = 0, we recover 
the pair contact process (PCP) introduced by Jensen [12| . 
In that case, any particle configuration without nearest- 
neighbour pairs is absorbing. The number of absorbing 
states grows exponentially with the chain length L and 
it is given by the Fibonacci number (see appendix A). 
The steady-state properties of the d = phase transi- 
tion between the active and inactive phases are described 
by the directed percolation (DP) universality class [12[. 
However, the dynamical properties of this transition are 
more subtle and still under active investigation p3|.p4[. 
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As we shall see, the presence of single-particle diffusion 
changes the universality class of the transition between 
the active and inactive phases: for any finite values of d 
it does not fall in the DP universality class anymore. 

Long ago, Janssen and Grassberger (25) conjectured 
that a model with a continuous phase transition from 
a fluctuating active phase into a phase with a single ab- 
sorbing state and without additional symmetries is in the 
DP universality class. While it is widely believed that in 
the presence of local symmetries in the reaction rates 
there should be a different universality class (normally 
the PC if the particle number is locally conserved mod- 
ulo 2), Park and Park J26j provided a counterexample 
which shows that even in the presence of local symmetries 
the steady-state transition can be in the DP universality 
class. This already indicates that the universality clas- 
sification might be more subtle than previously thought. 
The exponents we want to compare with are (see fj],§]) 

DP: 9 = i/y/i/x ^ 1.5806... , P/u± ~ 0.2520 .. . 
PC: = z/||/i/ x ^ 1.749... , 0/v x ^ 0.499 .. . (11) 

and are defined as usual from the density n(p) ~ (p—Pc)^ 3 
and the spatial and temporal correlation lengths ~ 
(p — p c )~ u±A] , see also section IV. 

The stochastic time evolution of the system is deter- 
mined by the master equation, cast into the form 

^ = -H|P(.» (12) 

where \P(t)) a state vector and H is referred to as "quan- 
tum" Hamiltonian. For a chain with L sites, H is a 
stochastic 2 L x 2 L matrix with elements 

(<t\H\t) = -w(t -> a) , (a\H\a) = J2w(a ^ t) 

where \a), |r) are the state vectors of the particle config- 
urations a, t and w are the transition rates. 

Since H is non-hermitian, it has distinct left and right 
eigenvectors. We will use the notation |0 r ), |l r ), . . . , \n r ) 
((0j|, (L|, . . . , (ni\) for the right (left) eigenvectors of H 
corresponding to energy levels Eq, E%, . . . , E n ordered 
according to diE < Si-Ei < . . -$tE n , where 3? denotes 
the real part. One has (ni\m r ) = if E n ^ E rn and we 
normalize the states in such a way that (ni\n r ) = 1. 

Steady states are right eigenvectors of H with zero 
eigenvalue. The two absorbing configurations mentioned 
above are steady states and given by 

|0 r ) := I000...0) (13) 
|l r ) := | A00 . . . 0> + |0A0 ... 0) + ... + |000 ... A) (14) 

Thus, H has two zero eigenvalues Eq — E\ =0, while 
r := ^{Ei is the inverse relaxation time towards the 
steady state. Furthermore, since H is stochastic, 



(0*h=EH ( 15 ) 

is a left eigenvector of H with zero eigenvalue. 

The calculation of the eigenvalues and eigenvectors of 
the stochastic Hamiltonian H, from which we will derive 
the critical properties of the model, is performed by the 
dmrg algorithm (l4|,[ll| adapted to non-hermitian matri- 
ces, as described in detail in In section III we present 
a further improvement to deal with systems with degen- 
erate ground states. 

III. SHIFT OF THE TRIVIAL GROUND STATE 

A severe numerical problem is generated by the fact 
that the physically relevant eigenstate, the first excited 
state, is only the third state in the spectrum due to the 
double degeneracy of the ground state. Asymmetric diag- 
onalization algorithms for large sparse matrices are much 
less robust than their counterparts for symmetric matri- 
ces, making a precise determination of eigenvalues not at 
the extreme ends of the spectrum much more demand- 
ing. In particular, the precision of the eigenstates suffers. 
This in turn leads to a lower quality truncation of the 
state space in the dmrg algorithm and subsequently to 
further deterioration of the results for longer chains. In 
the present case, the DMRG becomes unstable for rather 
short chain lengths (L ss 20). 

We alleviate this problem successfully by projecting 
out one of the two ground states, for which both the 
left and right eigenstate are known, as given in eqs. 
(HiDi)- All non-degenerate eigenstates of the asym- 
metric Hamiltonian H obey a biorthogonality condi- 
tion. This leaves us with two options to eliminate one 
ground state: (i) Each diagonalization algorithm for 
large matrices (Arnoldi or Lanczos) generates a sequence 
of (bi)orthogonal trial vectors, such that the requested 
eigenstate can be written as a linear combination of 
those. During the generation, one may enforce orthogo- 
nality of all these trial vectors with respect to the ground 
state. This option does not exist if one uses a black box 
routine, (ii) One may directly modify the Hamiltonian 
shifting the ground state to an arbitrarily high energy, 
i.e. working with 

H'(A) :=JJ + A|0 r )(0i| (16) 

where A is a positive number larger than the energy gap 
of H. This new Hamiltonian is not stochastic, however it 
has the same spectrum as H apart from one of the ground 
states which is shifted to an energy level A. The gap can 
now be obtained from the first excited state of H'(A). 
This procedure can also be implemented for a stochastic 
Hamiltonian with a single ground state. The gap is then 
obtained from the ground state energy of H' (A). This 
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option can always be used independently of the diago- 
nalization method employed, see Q in relation with the 
power method. 

To implement both options, one has to write down 
the left and right ground states in the transformed block 
bases generated by the DMRG . 

Let us denote for a block of length L by |0l) the state 
with no particles and by (t£ | = J2a L ( a L I the sum over all 
states of a block (i.e. in the complete, unreduced Hilbcrt 
space of the block). In the reduced block basis 
produced by the DMRG we have approximately |0i) = 
J2 mL ( m L\^L)\m L ) and (Ti| = J2 mL ( T L\m L )(m L \. The 
right and left ground eigenstates then read 

|0 r ) 2 i +2 = |0i) <8> |0) <8> |0) ® |0i) = 

(mi|0i)(mi|0i)(|mi)®|0)®|0)®K)), (17) 

rriL ,m' L 

(Oihi+2 = ® (s\ (8 (s'| ® (t l | = 

s,s' 

(TL\m L ){T L \m' L )({m L \®{s\®{s'\®{m' L \), (18) 

rriL ,s,s' 

where s and s' run over single site states. 

At the beginning of the DMRG application, the matrix 
elements (tol|0l) and (Ti|mi) are trivially constructed 
for a sufficiently small L such that the Hilbert space has 
less then m states, the number of states kept. 

Using (t l+1 \ = £ s (-n|®( s l and |0 L+ i) = |0i)®|0) and 
inserting one in (t l+ i| = E mL+1 ( T £+il m £+i)( TO i+il 
and |0i+i) = S mi+1 (mi+i|0i + i)|mi + i) one finds as 
recursive relation 

(mi+i|0i+i) =^(m L+1 |m L 0)(m L |0 L ), (19) 

rriL 

(T L+1 \m L+1 ) = ^2{m L s\m L+1 ){T L \m L ). (20) 

The matrix elements (mi+i|mis) are from the incom- 
plete basis transformation |mi) ® \s) — > |mi+i). 

Numerical implementation of both approaches reveals 
that the second approach is numerically very stable, 
while the first one is not, at least not if one carries 
out an unsophisticated Gram-Schmidt orthogonalisation. 
We suppose that this is due to the fact that global or- 
thogonality of the trial vectors is numerically not ex- 
actly conserved by the diagonalization algorithms, while 
Gram-Schmidt assumes this when a new trial vector is 
added and made orthogonal to the ground state. To find 
r = -JfE^, we have chosen as density matrix 

p = |tr{|0 r )(0 r | + |0,)(0,| + |2 r )(2 r | + |2,)<2j|} (21) 

where tr denotes the partial trace on the states of the 
left or right side of the chain. It is essential to target 
also the ground state |0 r ), (0;| projected out to maintain 



a good description of the ground state via (0) and JTs| ) 
after some DMRG steps. Otherwise, the Hamiltonian still 
contains a small perturbing contribution from the zero- 
energy ground state, which might destroy the stability of 
the procedure. 



IV. NUMERICAL RESULTS 
A. Analysis of the gap 

1. Finite-size scaling method 
We are interested in the lowest energy gap 



T = T{p, d- L) = E 2 



(22) 



whose finite-size scaling behaviour contains the desired 
information about the critical properties of the PCPD 
model (|I],|]j3]). While in principle the eigenvalues of the 
(non-symmetric) matrix H may have a non-zero imagi- 
nary part, we always found that Ei is real. 

Generically, we expect the following finite-size be- 
haviour of r 




if p < p c {d) 
if p=p c (d) 
if p > p c (d) 



(23) 



with 9 = v\\/v±_ and £j_ the spatial correlation length. 
At a continuous transition this diverges as ~ \p — 
p c {d)\~ v± , while in the time direction £|| ~ \p— Pc(d)\~ v ^ ■ 
The first line in ( |23| ) indicates that on a finite lat- 
tice, the relaxation towards the absorbing configurations 
is very slow. Indeed, for L — > oo the relaxation time 
t = r _1 becomes infinite and a state with a finite den- 
sity of particles becomes a steady state of the system. 
From the third line in (23) we see that in the entire in- 



active phase the scaling behaviour of the energy gap is 
algebraic, with an exponent equal to 2, as for the pair- 
annihilation model 2A — > 0. Therefore, in addition of 
the obvious double degeneracy of the ground state, we 
find that our model is gapless in the large-system limit, 
linii^oo T(p, d; L) = 0, for all values of p and d (in most 
systems the gap is finite in the inactive phase). 

The different phases and the critical point can be best 
identified from the analysis of the quantity 



Y L (p,d) := 



lnrr(p,rf;£ + l)/r(p,d;£-l)] 
In [(£ + !)/(£-!)] 



(24) 



While usually, the critical point is found by looking for 
intersections of two curves Yj, and Yjy for two different 
lattice sizes L, L', it turns out that because of the scaling 
( p3| ) , the critical point is found from the maximum of Yi 
as a function of p for fixed d and L. From (23) one then 
has for the large-L behaviour of Yi (p, d) 
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( -L/U ; if p < p c (d) 
Y L (p,d)~ \ -9 ; if p= Pc (d) (25) 
[ -2 ; if p>p c (d) 

Therefore, the location of the maximum value of Y^(p, d) 
yields a sequence of estimates p c (d; L) and the critical 
point can be found from extrapolating these sequences 
to L — ► oo. We point out, however, that the value of 
Yi, at its maximum cannot be used to infer the value of 
the exponent 9. Only after p c (d) is determined for the 
L — > oo lattice, we can use (|25|) to find 9. The technical- 
ities of the method are discussed in appendix B [^7j . 

2. Inactive phase 

Figure |l| shows a plot of Yl (p, d) as function of the pa- 
rameter p and for d = 0.5 for L = 9, 11, 13, 15 and 17 
(in the inset we show Yl for d = 0.2). Calculations were 
performed up to L = 30, but for the largest sizes only in 
the vicinity of the critical point. 



limit value. This check also confirms that the numerical 
values for the energy gap, as obtained from the dmrg 
calculation, are very accurate. 



TABLE I. Finite-size estimates — Yl of the dynamical ex- 
ponent 9 and their L — > oo extrapolation, obtained by the 
BST method, in the inactive phase. 



L 


d = 0.1, p = 0.6 


d = 0.5, p = 0.6 


d = Q.8, p = 0.5 


9 


1.979223431 


1.7110633 


1.226133000 


11 


1.986196898 


1.7567773 


1.267952170 


13 


1.990160149 


1.7893116 


1.306670859 


15 


1.992629239 


1.8139503 


1.342357293 


17 


1.994271913 


1.8333674 


1.375169517 


19 


1.995420176 


1.8491053 


1.405319671 


21 


1.996254444 


1.8621348 


1.433037552 


23 


1.996879687 


1.8731052 


1.458550029 


25 


1.997360512 


1.8905590 


1.482070754 


27 


1.995400349 


1.8976147 


1.503795737 


29 






1.523902146 


oo 


2.0000(1) 


2.000(1) 


1.99(2) 



-1.4 



>h J -1.6 



-1.8 



° I 9 

■ L-1 1 
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L=15 
< i L=17 


-1.8 

X -2 


d = - 2 


8 N> 0.1 0.15 0.2 
d = °- 5 Inactive 



0.2 0.4 0.6 0.8 1 

P 

FIG. 1. Plot of Yl(p, d) as function of p for various lattice 
sizes and d — 0.5. The inset shows the same quantity for 
d = 0.2. 

In the active phase (small p) the quantity Yl (p, d) 
scales linearly with L as predicted in (|29). For larger 
p the reaction 2A — > dominates and the system is in 
the inactive phase, for which we should have asymptot- 
ically Yl(p, d) ~ —2. As the system size L is increased, 
the curves approach the expected value —2 in the whole 
inactive phase. Table | shows the values of — Yl(p, d) cal- 
culated for three different diffusion constants and in the 
inactive phase. We extrapolated the finite-lattice data, 
using the BST algorithm [^8| and found excellent agree- 
ment with the expected value of 9 = 2, see eq. (p5|). The 
convergence is similar to the examples studied in al- 
though the raw data can be quite far from their L — > oo 



3. Active-inactive transition line 

Table || collects our results for the critical point p c (d) 
and the exponents 9 and /3/fj_, for several values of d. 



TABLE II. Critical parameters p c (d), and fl/vx along 
the active-inactive transition line for several values of d. 



d 


0.10 0.15 0.20 0.35 0.50 0.80 


Pc 

9 

Plvx 


0.111(2) 0.116(2) 0.121(3) 0.138(1) 0.154(1) 0.205(3) 
1.87(3) 1.84(3) 1.83(3) 1.72(3) 1.70(3) 1.60(5) 
0.50(3) 0.49(3) 0.49(3) 0.47(3) 0.48(3) 0.51(3) 



First, estimates for p c (d) in the L — > oo limit were 
obtained from a linear fit in l/L of the finite-size data 
p c {L), i.e. the abscissas of the maxima of Yl (p, d) . For 
comparison, we recall that p c {0) = 0.077090(5) @. 

Next, we estimate the value of the exponent 9 from 
the extrapolation of 9l = —YL{p c {d),d) to the thermo- 
dynamic limit L — > oo. Figure |^ shows a plot of 9l as 
a function of l/L for various values of d. We notice a 
different finite-size scaling behaviour of 9l for d < 0.20 
and for d = 0.50. In the former case 9l varies quite 
strongly with L, starting from a value above 2 and going 
towards values below 2. For d w 0.5, 9l shows little in- 
dependence. The extrapolated values of 9 (coming from 
a cubic fit in l/L) are shown in Table O. Error bars are 
due to the uncertainty in the determination of the critical 
point location. The final estimates for 9 appear to vary 
with d. It is not yet clear whether this variation is real 
or the consequence of an unresolved finite-size correction 
term. However, the results are clearly inconsistent with 
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a transition of DP type, for which 6 w 1.58. 



I d = 0.10,p = 0.111 
I d = 0.15,p = 0.116 
■ d = 0.20, p = 0.121 
' d = 0.35, p = 0.138 
d = 0.50, p = 0.154 
l d = 0.80, p = 0.205 
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FIG. 2. Plot of 6l as a function of 1/L for various values of 
d. The two horizontal arrows point to the values of 9 expected 
for the DP (9 » 1.58) and PC (9 « 1.75) universality classes. 



B. Steady state particle density 

For a finite lattice and for all values of p and d the sta- 
tionary states, given in ( |l3| ) and (|l4|), contain either no 
or just one particle. Therefore the steady-state particle- 
density vanishes in the large L limit. 

To obtain a steady state with a finite density of par- 
ticles we added a particle creation process at the two 
boundary sites M 



A with rate p , 



(26) 



We are interested in the (steady state) density profile, for 
a chain of length L 



n L (l) = (Q/|n(Z)|0 r ) 



(27) 



where I — 1,2, ... L labels the position along the chain, 
n(l) is the particle number operator at site I and |0 r ) 
and (0;| are the left and right ground states of the non- 
symmetric operator H, with the terms coming from j2^ ) 
added. 

The boundary reaction term removes the ground state 
degeneracy of H encountered earlier. Therefore, dmrg 
calculations are easier to perform when p' ^ 0. In this 
case it is not necessary to follow the shift strategy in- 
troduced in Sec. III. In fact the DMRG calculations are 
stable up to chains of lengths L ss 50 — 60, i.e. almost 
of the double size with respect of the lengths we could 
reach in the study of the energy gap. 

In general, the time-dependent particle density n — 
n(t,p;l,L) at site I for a lattice of size L and in the 



vicinity of the critical point, should satisfy the scaling 
form 

n(t, p- I, L) = F (t/£|| , L/U, l/L) 

"tl^Gitf&L/^l/L) (28) 

where for simplicity the dependence on d and p' is sup- 
pressed and F and G are scaling functions. The expo- 
nents have their usual meaning In particular, the 
steady-state density profile, at the critical point, should 
satisfy 



n L (l) := n(^, Pc ;l,L) = L~^f(l/L) 
with some scaling function f(z). 

30 



20 



(29) 



10 



0.2 0.4 0.6 0.8 

1/L 



FIG. 3. Scaled particle density Lnr,(l) as function of the 
scaled variable l/L in the inactive phase for d = 0.5, p = 0.8 
and injection rate p' = 0.3. Data for different system sizes 
(L — 20, 24, 28, . . . 48) collapse nicely onto a single curve. 

First, we consider the inactive phase. For our model, 
the average particle density decays algebraically for large 
times as n(t) ~ i~ 1/2 . From @, we identify P/u\\ = 1/2. 
Thus f3/v± — 0f3/v\\ — 1 since the anisotropy exponent 
9 = V\\/i>x = 2 in the inactive phase, see table |. 

The scaling (29) is confirmed in Fig. [| obtained for 
chains up to L = 48, which shows the scaled particle den- 
sity Ltil{1) as a function of l/L. The ratio f3/v± = 1 for 
the whole inactive phase has also been confirmed with 
good accuracy from BST extrapolations. 

From now on we focus on the transition between the 
active and inactive phases. We concentrate on the central 
particle density 



n{L) :=n L (L/2) 



(30) 



For this quantity we expect that for L large, n(L) ~ 
no + 0(e~ L ^ ± ) > in the active phase where at the 
critical point 



n(L) - L-VI^ 



(31) 
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while, as shown above, n(L) ~ L^ 1 in the entire inactive 
phase. 

Again, in analogy with ( |24| ) we form the logarithmic 
derivative 



\n[n(L + l)/n(L-l)} 
' ln[(L + l)/(L-l)] 



(32) 



According to eq. (|3l|), p(L) is expected to converge to 
the ratio for L — -> oo at the critical point p = p c (d) 

as listed in Table ||. A plot of p(L) for various values 
of the diffusion constant is shown in Fig. [l| We notice 
that at small L and d the effective exponent p{L) starts 
from values close to DP, while it clearly deviates from it 
for increasing L. A cubic fit in 1/L yields the extrapo- 
lated exponent ratio (3/v± as a function of d, see Table 0. 
This ratio appears to be rather constant with d. As in the 
case of the calculation of 9, we notice that the exponents 
are clearly inconsistent with a DP transition. However, 
the extrapolated values are close to the exponent value 
expected for a transition in the PC class. 
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FIG. 4. The effective exponent p(L) as function of 1/L for 
various values of d and up to L — 51. The horizontal arrows 
point to the values expected for DP and for PC universality 
classes. 



V. DISCUSSION 

We collect the results of this study of the Pair Contact 
Process with Diffusion. In Fig. ||(a) we show the steady- 
state phase diagram of the PCPD model. For small 
enough p, there is an active phase with a non-vanishing 
steady-state particle density and which goes over into an 
inactive state at some critical point p c (d). The critical 
line separating the active from the inactive phases termi- 
nates for d -> at the DP point p c {0) = 0.077090(5) [[L3j. 
For d — > 1 the critical line terminates at the MF point 
p = 1/3, as predicted from the mean field equations. 
This should have been expected, since it is well-known 



that when diffusion dominates over all other reactions 
the critical behaviour becomes of mean field type, even 
in one dimension 

In the PCPD, the entire inactive phase is critical and 
is expected to be in the universality class of diffusion- 
annihilation We have found the exponents 9 = 2 
and /3/v± = 1, confirming this expectation. 

We have investigated the properties of the transition 
line between the active and inactive phases. Fig. ||(b) 
and (c) show the numerical estimates of the exponents 9 
and f3 /v±, respectively, see also Table ||. For compari- 
son, the values ([ll]) of these exponents for the DP and PC 
universality classes are shown as horizontal lines. Clearly, 
our results are incompatible with a transition in the DP 
class. That means that single particle diffusion is a rel- 
evant perturbation of the pair contact process. While 
(3/v± is quite constant and consistent with the PC value 
within error bars, 9 seems to vary continuously with d. 
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FIG. 5. (a) Steady-state phase diagram of the model 
2A — » 3A, 2A — > as determined by dmrg techniques. The 
continuous line separates the active from the inactive phase. 
Extrapolated estimates of the exponent 6 are shown in (b) 
and of the exponent /3/u± in (c) as a function of d. 

There is no clear evidence for two distinct universality 
classes along the critical line, as predicted by pair mean 
field theory (see section II), but our information on the 
transition line for d close to unity remains incomplete 
(see also below). At present, neither of the two possi- 
bilities can be ruled out. It is an open question whether 
the 2D PCPD model phase diagram will be in qualitative 
agreement with pair mean field theory. 

Turning to the apparent variation of 9 with d, note that 
decreases systematically with d, starting from 9 « 2 in 
the d — > limit. To explain this, consider the case when 
d is small and the particle density is low, i.e. the system 
is close to criticality. Then the particles will diffuse inde- 
pendently until they meet and react. The lower the diffu- 
sion rate d, the larger the time interval tdis for which free 
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diffusion can be observed. Associated to this is a length 
scale -Rdiff, viz. tdiff ~ Rdiff- For times short compared to 
tdie, the system is governed by diffusion only. Therefore, 
for system sizes L <§C Rdiff , the critical exponents will be 
effectively those of free diffusion and 8 ~ 8dm = 2. The 
true value for 8 will be seen only if L ^> Rdie- These 
two regimes are separated by a crossover region and it 
appears plausible that the apparent variation of our es- 
timates of 8 with d might be the consequence of such 
a crossover. However, this crossover phenomenon does 
not show up in the calculation of f3/v±. That is to be 
expected since by injecting particles at the boundaries 
we induce a finite density of particles also at the critical 
point and tdie will no longer increase beyond any bounds 
as d — > and the above heuristic argument no longer 
applies. 

While the exponent (3/v±_ is not far from the one of the 
PC class it is difficult to extract reliable information form 
8. We have no reason to believe that exponents should 
be continuously varying as function of the diffusion con- 
stant. A closer inspection of the finite-size behaviour (see 
Fig. ||) reveals that the L — > oo asymptotic value 8 is ap- 
proached from above for d < 0.35, while it is approached 
from below for d > 0.5. As finite-size effects are seen 
to be weak in the interval 0.35 < d < 0.5, we presume 
that the most reliable estimates for 8 might be those 
obtained in this range of diffusion constants. Therefore 
most likely 8 « 1.7, which is not far from the PC value 
either. We also recall that the relaxation in the inactive 
phase is algebraic: All known models |26|,|36|-[39| in the 
PC class are characterized by an algebraic relaxation in 
the inactive phase with exponents 8 = 2 and (3/v±_ = 1, 
as for the PCPD model. Therefore, it might be tempt- 
ing to speculate that the active- inactive transition of the 
PCPD model were in the PC universality class. 

On the other side all known models in the PC univer- 
sality class are characterized by some conservation laws 
either on the parity of the number of particles or by an 
exact symmetry between the absorbing states ||. Such 
local consevation laws are however lacking for the PCPD 
model, which suggests that the PCPD model should not 
be in the PC class. It should be also stressed that the 
unambigous identification of a steady-state universality 
class requires the determination of four independent ex- 
ponents, see J^||, while our techniques provided values 
for only two. 

In summary, we have used the dmrg method to find 
the steady-state phase diagram of the PCPD model. Sin- 
gle particle diffusion is a relevant perturbation for the 
system since the model does not show a DP behaviour as 
in the limit of vanishing diffusion constant. 

Some exponent values along the active-inactive tran- 
sition line are surprisingly close to those of the PC uni- 
versality class. However, it is conceivable that the near 
coincidence of the exponents with those of the PC class 
might be accidental. In that case the transition(s) would 



belong to a universality class distinct form both DP and 
PC. It is not yet clear whether in the ID PCPD there are 
two distinct transitions, as suggested by pair mean field 
theory, or merely a single one. All in all, complementary 
studies would be needed to fully understand the remark- 
ably subtle behaviour of this so simple-looking model. 

Note added: After this paper was submitted, Hinrich- 



scn 1 40 performed a Monte Carlo study of the PCPD for 



d = 0.1. The time-dependent density n(t) ~ t is char- 
acterized by the exponent 6 = (3/{v±8). He finds 6 = 
0.25(2) and 8 = 1.83(5). This agrees with our Tabic ||. 
After this paper was accepted, we also received a paper 
by Odor jtl| which studies the PCPD through Monte 
Carlo simulations and the coherent anomaly method. In 
particular, for 0.05 < d < 0.2, Odor finds S w 0.27 and for 
d = 0.5 and d = 0.9, S w 0.2. The first result compares 
well with our results from Table |J, which gives S w 0.27 
and agrees with the result of Hinrichsen POJ. However, 
Odor also finds « 0.58 for 0.05 < d < 0.2, consistent 
with the upper bound [3 < 0.67 reported_in_[|40| . That is 
far away from the PC value [3 PC — 0.93 | 
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APPENDIX A: NUMBER OF ABSORBING 
CONFIGURATIONS IN THE PCP 

For d = 0, the model reduces to the pair contact pro- 
cess (PCP) p2[ . In this case, all configurations of the 
type | . . . 0A0 . . . 0A0 . . .), without nearest neighbour par- 
ticles are absorbing and stationary states. We find the 
number N(L) of absorbing states in the PCP for a chain 
of length L with free (periodic) boundary conditions. 

For free boundary conditions, one has the recursion 



N(L) = N(L - 1) + N(L - 2). 



(33) 



To see this, concentrate on the leftmost site. If it is oc- 
cupied, its neighbour must be empty for the state to be 
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absorbing and there remains an open chain of L — 2 sites 
to be considered. If it is empty, one considers the open 
chain of the remaining L — l sites. The initial conditions 
for the problem are N(l) = 2 and N(2) = 3. Therefore 
N(L) = F L+1 is the (L + l) th Fibonacci number. 
This can also be seen from the generating function 



£ 

k=Q 



N(s)^Y j N(k)s k 



1 + s 



which satisfies because of ( |33] ) 

- _ N(0)(l- s ) + N(l)s 

[S > ~ 1 - a - s 2 ~ 1 - s - s 2 

For the inverse transformation one can use 

i r , A , , _ L _! 



iV(L) = — * ds N(s) s 

ZTTl 



(34) 



(35) 



(36) 



where the integral is taken in a closed circle centered in 
the origin of the complex s-plane and with radius smaller 
than the radius of convergence of the series (|34j) , so as to 
exclude contributions of other poles than that in s = 0. 
The transformation z = 1/s yields 



N(L) = 



1 

2tti 



1 + z L 
dz- -z L = 

z z — z — 1 



z L+2 _ z L+2 



(37) 



where z± = (1 ± v5)/2; z+ is the golden mean. 

For periodic boundary conditions, the number N peT (L) 
of absorbing states for a chain of L sites is 

N pcr (L) = N(L - 1) + N(L -Z) = z\ + z L _ (38) 

For L large the asymptotic behaviour is N(L) ~ 
N pcr (L) ~ z %. 

Another example of a kinetic model with exponentially 
many absorbing states is an adsorption-desorption model 
of fc-mers (fc > 3) on a ID lattice where the number 
of absorbing states is described by generalised Fibonacci 
numbers IH3I. 



APPENDIX B: FINITE-SIZE-SCALING IN 
SYSTEMS WITH A CRITICAL PHASE 

We discuss the finite-size scaling of the lowest gap 
T = Tl(j>) and how to localize the critical point. For 
simplicity, we suppress the dependence on d and write 
the scaling form 



r L (p) = L- e f((p- Pc )L 1 /^) 



(39) 



where / is assumed to be continously differentiable. For 
the gap, one expects the behaviour 



Tl(p) 



e -aL 

L- 2 
r 

1 oo 



if P < Pc 
ifp> pc, 
ifp> pc, 



case C 
case N 



(40) 



where <7, are constants independent of L. Here, case 
N refers to the "normal" case of non-critical phases on 
both sides of the transition at p = p c and case C alludes 
to the case of a critical phase on one side and we have 
already set the exponent equal to 2 in view eq. ( p3[ ) valid 
for our model. This implies for the scaling function f(z) 



exp(-A\z\"±) 

(fi-2)v x 



9u± 



z — * — oo 

z — * +oo, case C (41) 
z — * +oo, case N 



where A is a positive constant. Since f(z) is positive, it 
follows that in the case C with 9 < 2, f(z) must have 
a maximum at some finite value z max - For the case N, 
however, f(z) should increase monotonically with z. 
The estimator Yj, as defined in (E4j) then becomes 



Y L = -t 



ln[/(z + )//(z_)] 
ln[(£ + l)/(L-l)] 



(42) 



where z± — (p — p c )(L ± l) 1 ^ u± . Furthermore, writing 
g(z) — ln/(z), a straightforward calculation gives 



lim • 



[g'(z)-zg"(z)] 



dY L _ L 1 /^ 
dp is ± 



(43) 



oo 



in the finite-size scaling limit p — > p c and L — > oo simul- 
taneously such that z — (p ~ p c )L 1 / 1/± is kept fixed. For 
the case C with 9 < 2 and v± < 2, there is some finite z* 
such that dY L /dp\ z=z , = 0. Then 



Y L {z*) = -t 



+ —z*g'(z*) 



(44) 



If we choose z = z*, we have a sequence of values of p L 
converging towards p c according to p^ ~ p c + z*L~ 1 ^ u± . 
On the other hand, p^ can be found by determining the 
maximum oJYl as a junction ofp, since limdYL/dp| z » = 
0. This is the desired result. However, because of (|44|) 
Yl(z*) does not readily yield an estimate for the expo- 
nent 9, since there is no guarantee that g'{z*) should 
vanish, (or in other words that z* = z max ). The gener- 
alization to other observables with scaling analogous to 
( p0| ) is immediate. 

For the case N, however, this technique does not apply, 
since in general f'(z) > for all values of z. 

Finally, we recall that the leading finite-size correction 
terms determine whether or not the curves Y^(p) and 
Yl'(j>) will intersect. Consider the extended scaling form 
T L (p) = L- f(z) [1 + L~ u A(z)] , where u > is the lead- 
ing correction exponent. If uj < 2, we find 



Yr = -6 



Z9,{z Kl-»A(z)(-u, + ± Z 4^) (45) 



v±_ A{z) 
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up to terms of order 0(L~ 2 , L^ 1 ^^). Now, the curves 
Yt, and Yl> intersect if there is some z; nt such that the 
scaling function of the leading correction term in ( |45| ) 
vanishes. But that term depends only on the correction 
amplitude A(z) and is independent of f{z). 



Unite Mixte de Recherche CNRS No 7556 
[1] V. Privman (Ed.) Nonequilibrium Statistical Mechanics 
in One Dimension, Cambridge University Press (Cam- 
bridge 1996) 

[2] J. Marro and R. Dickman, Nonequilibrium Phase Tran- 
sitions in Lattice Models, Cambridge University Press 
(Cambridge 1999) 

[3] H. Hinrichsen, Adv. Phys. 49, 1 (2000) 

[4] J. Cardy and U. C. Tauber, Phys. Rev. Lett. 77, 4780 
(1996) and J. Stat. Phys 90, 1 (1998); D.C. Mattis and 
M.L. Glasser, Rev. Mod. Phys. 70, 979 (1998) 

[5] D. Kandel, E. Domany and B. Nienhuis, J. Phys. A: 
Math. Gen. 23, L755 (1990); F.C. Alcaraz, M. Droz, 
M. Henkel and V. Rittenberg, Ann. of Phys. 230, 250 
(1994); I. Peschel, U. Schultze and V. Rittenberg, Nucl. 
Phys. B430, 633 (1995); P.-A. Bares and M. Mobilia, 
Phys. Rev. Lett. 83, 5214 (1999) and 85, 893 (2000); 
CM. Schiitz, Integrable Stochastic Many-Body Systems, 
to appear in C. Domb and J.L. Lebowitz (Eds.), Phase 
Transitions and Critical Phenomena, Vol. 19, Academic 
Press (New York 2000). 

[6] B. Chopard and M. Droz, Cellular Automata Modelling 
of Physical Systems, Cambridge University Press (Cam- 
bridge 1998) 

[7] M. Henkel and H. Herrmann, J. Phys. A: Math. Gen. 23, 
3719 (1990); J.R.G. de Mendonca, Phys. Rev. E 60, 1329 
(1999). 

[8] J.R.G. de Mendonca, J. Phys. A Math. Gen. 32, L467 
(1999). 

[9] E. Carlon, M. Henkel and U. Schollwock, Eur. Phys. J. 

B12, 99 (1999). 
[10] J. Hooyberghs and C. Vanderzande, J. Phys. A: Math. 

Gen. 33, 907 (2000). 
[11] M.J. Howard and U.C. Tauber, J. Phys. A: Math. Gen. 

30, 7721 (1997). 
[12] I. Jensen, Phys. Rev. Lett. 70, 1465 (1993). 
[13] R. Dickman and J. Kamphorst Leal da Silva, Phys. Rev. 

E58, 4266 (1998). 
[14] S.R. White, Phys. Rev. Lett. 69, 2863 (1992); Phys. Rev. 

B 48, 10345 (1993). 
[15] I. Peschel, X. Wang, M. Kaulke and K. Hallberg (Eds.): 

Density Matrix Renormalization: A New Numerical 

Method in Physics, Springer, (Heidelberg 1999). 
[16] R.J. Bursill, T. Xiang and GA. Gehring, J. Phys. Cond. 

Mat. 8, L583 (1996); X. Wang and T. Xiang, Phys. Rev. 



B56, 5061 (1997); N. Shibata, J. Phys. Soc. Jpn. 66, 2221 

(1997); K. Maisinger and U. Schollwock, Phys. Rev. Lett. 

81, 445 (1998). 
[17] Y. Hieida, J. Phys. Soc. Jpn. 67, 369 (1998). 
[18] M. Kaulke and I. Peschel, Eur. Phys. J. B5, 727 (1998). 
[19] 0. Legeza, M. Kaulke and I. Peschel, Ann. der Physik 8, 

153 (1999) 

[20] J. Kondev and J. B. Marston, Nucl. Phys. B 497, 639 
(1997); J. B. Marston and Shan-Wen Tsai, Phys. Rev. 
Lett. 82, 4906 (1999); T. Senthil, J. B. Marston, and M. 
P. A Fisher, Phys. Rev. B 60, 4245 (1999); Shan-Wen 
Tsai and J. B. Marston, Ann. der Physik, 8, 261 (1999). 

[21] L. Peliti, J. Phys. A: Math. Gen. 19, L365 (1986). 

[22] B.P. Lee, J. Phys. A27, 2633 (1994). 

[23] I. Jensen and R. Dickman, Phys. Rev. E 48, 1710 (1993); 
J.F.F. Mendes, R. Dickman, M. Henkel and M.C. Mar- 
ques, J. Phys. A: Math. Gen. 27, 3019 (1994); M.A. 
Munoz, G. Grinstein, R. Dickman and R. Livi, Phys. 
Rev. Lett. 76, 451 (1996); P. Grassberger, H. Chate and 
G. Rousseau, Phys. R ev. E 55, 2488 (199 7). 



[24 
[25 

[26 
[27 

[28 



[29 

[30 

[31 
[32 

[33 
[34. 
[35 

[36 



[37 

[38 



[39 
[40 
[41 
[42 



R. Dickman, preprint cond-mat/9909347 



H.K. Janssen, Z. Phys. B42, 151 (1981); P. Grassberger, 
Z. Phys. B47, 365 (1982) 

H. Park and H. Park, Physica A221, 97 (1995). 

We thank H. Hinrichsen for a useful discussion on that 

point. 

R. Burlisch and J. Stoer, Numer. Math. 6, 413 (1964); 
M. Henkel and G. Schiitz, J. Phys. A: Math. Gen. 21, 
2617 (1988). 

R. Ziff, E. Gulari and Y. Barshad, Phys. Rev. Lett. 56, 
2553 (1986). 

R. Bidaux, N. Boccara and H. Chate, Phys. Rev. A 39, 
3094 (1989). 

R. Dickman and T. Tome, Phys. Rev. A 44, 4833 (1991). 
G. Odor, N. Boccara and G. Szabo, Phys. Rev. E 48, 
3168 (1993). 

N. Menyhard and G. Odor, J. Phys. A28, 4505 (1995). 

G. Odor and A. Szolnoki, Phys. Rev. E 53, 2231 (1996). 
K. Oerding, F. van Wijland, J.-P. Leroy and H.J. Hil- 
horst, J. Stat. Phys. 99, 1365 (2000). 

P. Grassberger, F. Krause and T. von der Twer, J. Phys. 
A: Math. Gen. 17, L105 (1984); P. Grassberger, J. Phys. 
A: Math. Gen. 22, L1103 (1989). 

N. Menyhard, J. Phys. A: Math. Gen. 27, 6139 (1994). 
N. Inui and A. Yu. Tretyakov, Phys. Rev. Lett. 80, 5148 
(1998); M.C. Marques and J.F.F. Mendes, Eur. Phys. J. 
B12, 123 (1999). 

H. Hinrichsen, Phys. Rev. E 55, 2 19 (1997). 
H. Hinrichsen, 



cond-mat/0001177 



G. Odor, Phys. Rev. E62, R3027 (2000). 

M. Barma, M.D. Grynberg and R.B. Stinchcombe, Phys. 

Rev. Lett. 70, 1033 (1993). 



10 



